;+
; NAME:
;       HISA
;
; PURPOSE:
;       Generate a spectrum for sources in the IGPS overlap by
;       integrating over the source and subtracting off a boundary spectrum.
;
; CATEGORY:
;       distance-omnibus
;
; CALLING SEQUENCE:
;       spectrum = HISA( struct [,V_STD=v_std][,SPECTRUM=spectrum]
;                       [,BGPS_ML=bgps_ml][,HI_DIRS=hi_dirs][,BGPS_DIR=bgps_dir]
;                       [,MAP_DIR=map_dir][,OBJ_DIR=obj_dir][,/PLOT]
;                       [,ONSPEC=onspec][,BDRSPEC=bdrspec][,OBJMASK=objmask]
;                       [,BORDER=border][,X_ORIG=x_orig][,Y_ORIG=y_orig]
;                       [,X_BDR=x_bdr][,Y_BDR=y_bdr][,BHDR=bhdr][,R=r] )
;
; INPUTS:
;       STRUCT   -- BGPS source structure (see READ_BGPS_CSV), which
;                   includes source longitude and latitude
;       BGPS_ML  -- BGPS map locations IDL save file.  Defaults to:
;                   ./local/bgps_v102_map_locations.sav 
;
; OPTIONAL INPUTS:
;       HI_DIRS  -- Array of directories containing IGPS spectral line cubes.
;       BGPS_DIR -- Directory containing the BGPS data
;       MAP_DIR  -- Directory containing the BGPS IMAGE files
;       OBJ_DIR  -- Directory containing the BGPS object files
;       R        -- Rind radius, in BGPS pixels, for averaging
;                   the "off spectrum".  [Default: 12]
;
; KEYWORD PARAMETERS:
;       PLOT     -- Plot the spectrum to the screen
;
; OUTPUTS:
;       SPECTRUM -- Integrated masked spectrum for the BGPS source
;
; OPTIONAL OUTPUTS:
;       V_STD    -- Velocity vector associated with the SPECTRUM
;       SPECTRUM -- Integrated masked spectrum for the BGPS source --
;                   may be different than RETURN'd spectrum if
;                   routine hits an early exit.
;       ONSPEC   -- "On-Source" HI spectrum weighted by BGPS emission.
;       BDRSPEC  -- "Border" HI spectrum, averaged over admissible
;                   border pixels
;       OBJMASK  -- BGPS catalog object mask -- pixels from BGPS
;                   assigned to this object.
;       BORDER   -- "Rind" around the BGPS source, excluding pixels
;                   assigned to other catalog objects.
;       X_ORIG   -- X-axis pixel numbers for the BGPS source
;       Y_ORIG   -- Y-axis pixel numbers for the BGPS source
;       X_BDR    -- X-axis pixel numbers for the border region
;       Y_BDR    -- Y-axis pixel numbers for the border region
;       BHDR     -- BGPS image header
;
; COMMON BLOCKS:
;       HISA_BLOCK
;
; NOTES:
;    HI (neutral hydrogen) self-absorption can be used as a near/far
;    indicator If an absorption signature is seen at the velocity of
;    the source, it is likely to be at the near distance.  HI
;    self-absorption occurs because the galaxy is (approximately)
;    filled with warm (~100K) HI at all spatial positions and
;    therefore velocities.  A cold molecular cloud has enough neutral
;    HI column to be optically thick
;    (e.g. http://adsabs.harvard.edu/abs/2009ApJ...690..706A as a
;    recent example, but this was shown in earlier papers....)
;
;    Method: 
;    1. Acquire VGPS HI spectrum at given velocity
;       a. Check for the presence of a continuum source.  If one is
;          present, drop the weight of the HISA probability by a lot
;    2. Fit "baseline" at +/- 10 km/s.  
;       a. Set a weight based on how flat this baseline is
;    3. Measure the depth of absorption at specified velocity
;    4. Make a step-function similar to that of IRDCs with step depth
;       related to absorption well-depth 
;
; MODIFICATION HISTORY:
;    Created:  08/10/10, TPEB -- Initial Version, heavily modified
;                                from grsmatch.pro
;    Modified: 04/27/11, TPEB -- Code and documentation cleanup.
;    Modified: 04/29/11, TPEB -- Additional code and documentation
;                                cleanup, plus NO LONGER returning the
;                                square of the (on-off) spectrum, but
;                                simply ABS(on-off) (still enforcing
;                                ABSORPTION ONLY spectral features, as
;                                before).
;    Modified: 06/15/11, TPEB -- Replaced DELVARX with UNDEFINE as the
;                                means for clearing the data cubes
;                                from memory before loading the next
;                                one.  Fixed memory leak & improved
;                                reliability.
;    Modified: 06/27/11, TPEB -- Adjusted the masking for large +v in
;                                the VGPS and CGPS surveys (there are
;                                artifacts in the cubes at
;                                non-physical velocities).
;
;-

FUNCTION HISA, s, V_STD = v_std, SPECTRUM = spectrum, BGPS_ML = bgps_ml, $
               HI_DIRS = hi_dirs, BGPS_DIR = bgps_dir,$
               OBJ_DIR = obj_dir, MAP_DIR = map_dir, PLOT=plot,$
               onspec=onspec, bdrspec=bdrspec, objmask=objmask,$
               border=border, x_orig=x_orig, y_orig=y_orig,$
               x_bdr=x_bdr, y_bdr=y_bdr, bhdr=bhdr, R=r
  
  COMMON HISA_BLOCK
  
  ;; If the COMMON HISA_BLOCK is not defined by the calling routine,
  ;; then do things the long way 
  IF ~ exist(lastfn) THEN lastfn = ''
  
  ;; Set directory structure, if not specified in input
  IF n_elements(hi_dirs) EQ 0 THEN $
     hi_dirs = ['./local/sgps/','./local/vgps/','./local/cgps/']
  IF n_elements(bgps_dir) EQ 0 THEN bgps_dir = './local/bgps/'
  IF n_elements(obj_dir) EQ 0 THEN obj_dir = bgps_dir+'label/'
  IF n_elements(map_dir) EQ 0 THEN map_dir = bgps_dir+'maps/'
  
  ;; Define general parameters
  defsysv, '!MW', exists = exists
  IF NOT exists THEN galactic_params 
  V0 = !mw.v0
  
  ;; Define the velocity vector for use with the smoothed spectra
  v_std     = findgen(!MW.NVBINS)*!MW.DELTAV + !MW.VSTART
  spectrum  = fltarr(!MW.NVBINS)
  null_spec = spectrum
  onspec    = null_spec
  bdrspec   = null_spec
  
  ;; Read in Map Locations save file
  IF n_elements(bgps_ml) EQ 0 THEN $ 
     restore,file='./local/bgps_v102_map_locations.sav' ELSE $
        restore,file=bgps_ml
  
  ;; Select which BGPS object we're working with here
  hit = where(s.cnum eq bgps.cnum,ct)
  IF ct NE 1 THEN BEGIN
     message,'No match between distance structure and bgps SAVE file...',/inf
     RETURN, null_spec
  ENDIF
  bgps = bgps[hit]
  
  catroots = bgps.filename
  
  ;; Initialize structuring element for border mask
  ;; For GRS, r = 4 (1'), but HISA is a broader emission
  ;; structure, so take backgound as average over 3'.
  ;; r = 4      ; 1'  
  ;; r = 8      ; 2'  radius of this many arcmins
  
  IF n_elements(r) EQ 0 THEN r = 12
  
  elt = shift(dist(2*r+1,2*r+1),r,r) LE r
  
  
  ;; Check to see if this BGPS source lies within the IGPS coverage
  good = 0b
  IF (bgps.glon GE 4.53 AND bgps.glon LE 175) THEN good = 1b     ;; SGPS-CGPS
  IF (bgps.glon GE 350  AND bgps.glon LE 358) THEN good = 1b     ;; SGPS
  
  IF ~ good THEN BEGIN
     message,'BGPS source lies outside IGPS HI coverage...',/inf
     RETURN, null_spec
  ENDIF
  
  ;;=================================================================
  ;; Read in BGPS map and labelmap
  ;; data = BGPS map
  ;; obj  = Label map
  
  data = readfits(map_dir+catroots, bolo_hd, /SILENT)
  extast, bolo_hd, bolo_astrom
  sz = size(data)
  char = stregex(catroots, '_map')
  labelname = strmid(catroots, 0, char+1)+$
              'label'+strmid(catroots, char+1, 30)
  obj = readfits(obj_dir+labelname+'.gz', bhdr,/SILENT)
  
  
    ;; Determine where BOLOCAM points land in IGPS data
  objmask = obj EQ bgps.cnum
  
  ind1 = where(objmask, ct)
  IF ct EQ 0 THEN BEGIN
     message,'No match in BGPS label map...',/inf
     return, null_spec
  ENDIF
  
  ;; Assign weights & determine x,y coordinates of BGPS pixels
  wts = data[ind1]
  x_orig = ind1 mod sz[1]
  y_orig = ind1 / sz[1]
  
  ;; Define border region w/ elt and to exclude other BGPS sources
  border = (dilate(objmask, elt)-objmask)*(obj eq 0)
  ind2 = where(border,ctbdr)
  x_bdr = ind2 mod sz[1]
  y_bdr = ind2 / sz[1]
  
  ;; Get galactic coordinates for each pixelin BGPS source & border
  xy2ad, x_orig, y_orig, bolo_astrom, glon, glat
  xy2ad, x_bdr, y_bdr, bolo_astrom, glonbdr, glatbdr
  
  ;; Sort everything in order of increasing glon to minimize file reads
  wtorder = sort(glon)
  glon    = glon[wtorder]
  glat    = glat[wtorder]
  wts     = wts[wtorder]
  x_orig  = x_orig[wtorder]
  y_orig  = y_orig[wtorder]

  
  wt2order = sort(glonbdr)
  glonbdr  = glonbdr[wt2order]
  glatbdr  = glatbdr[wt2order]
  x_bdr    = x_bdr[wt2order]
  y_bdr    = y_bdr[wt2order]

  ;; Get filenames for HI cube and continuum map (& which survey)
  fn = get_hi_cube(glon, bgps.glat, hi_dirs, fnc, survey)
  fn2 = get_hi_cube(glonbdr, bgps.glat, hi_dirs, fn2c, survey2)
  
  survey  = mode(survey)
  survey2 = mode(survey2)
  
  ;; Build empty arrays for the smoothed spectra
  runspec = fltarr(!MW.NVBINS)
  bdrspec = fltarr(!MW.NVBINS)
  
  
  ;;====================================================================
  ;; Loop through all BGPS pixels assigned to this source
  FOR j = 0, n_elements(ind1)-1 DO BEGIN
     IF ~ file_test(fn[j],/READ) THEN CONTINUE
     
     ;; Pick HI filename and load a new one if necessary
     IF fn[j] NE lastfn THEN BEGIN
        UNDEFINE,hi_data
        print,fn[j]
        hi_data   = readfits(fn[j], hd)
        lastfn = fn[j]
        IF file_test(fnc[j],/READ) THEN BEGIN
           ;;print,fnc[j]
           cont_data = readfits(fnc[j], c_hd,/silent)
           have_cont = 1b
        ENDIF ELSE have_cont = 0b
     ENDIF
     
     ;; Get HI pixel coordinates for this BGPS pixel
     extast, hd, astrom
     rdhd, hd, s = h
     sz_hi   = size(hi_data)
     ad2xy, glon[j], glat[j], astrom, x, y
     
     IF have_cont THEN BEGIN
        extast, c_hd, castrom
        sz_cont = size(cont_data)
        ad2xy, glon[j], glat[j], castrom, cx, cy
     ENDIF
     
     ;; Look up spectrum -- Check first that (x,y) is within bounds of hi_data
     IF x LT 0 OR x GT sz_hi[1]-1 OR $
        y LT 0 OR y GT sz_hi[2]-1 THEN BEGIN
        wts[j] = 0
        CONTINUE
     ENDIF

     IF have_cont THEN BEGIN
        ;; Look at continuum -- if continuum emission T > limit, abort HISA
        IF cont_data[cx,cy] GE !MW.TA_CONT THEN BEGIN
           message,'21-cm continuum emission w/in BGPS source',/inf
           return,null_spec
        ENDIF 
     ENDIF
     
     spectrum = interpolate(hi_data,replicate(x,sz_hi[3]),$
                            replicate(y,sz_hi[3]),indgen(sz_hi[3]))
     
     ;; Interpolate to common velocity scale
     specout = interpol(spectrum, h.v, v_std)
     outside = where(v_std LT min(h.v) OR v_std GT max(h.v), nout)
     IF nout GT 0 THEN $
        specout[outside] = !values.f_nan
     runspec = runspec+specout*wts[j]
     
  ENDFOR   ;; End of loop for ON-SOURCE spectra extraction
  
  
  ;;====================================================================
  ;; Loop through all BORDER pixels
  FOR j = 0,n_elements(ind2)-1 DO BEGIN
     IF ~ file_test(fn2[j],/READ) THEN CONTINUE
     
     ;; Pick HI filename and load a new one if necessary
     IF fn2[j] NE lastfn THEN BEGIN
        UNDEFINE,hi_data
        print,fn2[j]
        hi_data   = readfits(fn2[j], hd)
        lastfn = fn2[j]
        IF file_test(fn2c[j],/READ) THEN BEGIN
           ;;print,fn2c[j]
           cont_data = readfits(fn2c[j], c_hd,/silent)
           have_cont = 1b
        ENDIF ELSE have_cont = 0b
     ENDIF
     
     ;; Get pixel coordinates
     extast, hd, astrom
     rdhd, hd, s = h
     sz_hi   = size(hi_data)
     ad2xy, glonbdr[j], glatbdr[j], astrom, x, y
     
     IF have_cont THEN BEGIN
        extast, c_hd, castrom
        sz_cont = size(cont_data)
        ad2xy, glonbdr[j], glatbdr[j], castrom, cx, cy
     ENDIF
        
     ;; Look up spectrum, with same check
     IF x LT 0 OR x GT sz_hi[1]-1 OR $
        y LT 0 OR y GT sz_hi[2]-1 THEN BEGIN
        ctbdr = ctbdr - 1
        CONTINUE
     ENDIF

     IF have_cont THEN BEGIN
        ;; Look at continuum -- if continuum emission exceeds T limit,
        ;;                      neglect this data point 
        IF cont_data[cx,cy] GE !MW.TA_CONT THEN BEGIN
           border[x_bdr[j],y_bdr[j]] = 0b
           ctbdr -= 1
           ;;print,cont_data[cx,cy]
           CONTINUE
        ENDIF
     ENDIF
     
     spectrum = interpolate(hi_data,replicate(x,sz_hi[3]),$
                            replicate(y,sz_hi[3]),indgen(sz_hi[3]))
     
     ;; Interpolate to common velocity scale
     specout = interpol(spectrum, h.v, v_std)
     outside = where(v_std LT min(h.v) OR v_std GT max(h.v), nout)
     IF nout GT 0 THEN $
        specout[outside] = !values.f_nan
     bdrspec = bdrspec+specout
     
  ENDFOR   ;; End of loop for OFF-SOURCE spectra extraction
  
  ;;==============================================================
  ;; Filter the spectra & create averaged spectra
  filter = savgol(32,32,0,6)
  
  specout = runspec/total(wts)
  spectrum = convol(specout,filter,/edge_trun,/nan) >0.
  
  bdrspec = bdrspec/ctbdr
  bdrspec = convol(bdrspec,filter,/edge_trun,/nan) >0.
  
  ;; Mask obviously spurious data (T_A > 300 K)
  spectrum *= (spectrum LE 300)
  bdrspec  *= (bdrspec LE 300)
  
  ;; Mask out the region at V_LSR >= V_MAX + 40 km/s.  There is
  ;; sometimes crazy noise in this part of the spectrum that affects
  ;; the HISA output spectrum.  For QUADRANT I, V_MAX = V_TAN, for
  ;; QUADRANT II, V_MAX = 0.  Additionally, this only seems necessary
  ;; for VGPS and CGPS data -- SGPS data does not seem to suffer from
  ;; these random fluctuations at large positive velocity.
  
  IF survey GT 1 THEN BEGIN            ;; l > 18 in QUADs I and II
     
     IF s.glon_peak LE 90. THEN BEGIN  ;; QUADRANT I
        
        ;; Get the TANGENT VELOCITY in the LSR frame (as opposed to the
        ;; physical radial velocity).
        d     = (dindgen(!MW.NBINS)*!MW.BINSIZE + !MW.BINSTART)
        vlsr  = vphys2vlsr(s.glon_peak, s.glat_peak, d)
        v_max = (INTERPOL( vlsr, d, !MW.R0*cos(s.glon_peak*!dtor), /LSQUAD))[0]
        
     ENDIF ELSE v_max = 0.d            ;; QUADRANT II
     
     vmask = v_std LE (v_max + 40.d)
  ENDIF ELSE BEGIN
     vmask = byte(v_std *0. + 1.)    ;; vmask is unity for SGPS
     v_max = 1.d3                    ;; Define as junk only if trying to plot
  ENDELSE                            ;; SGPS objects
  
  ;; Create the Final Spectrum     
  onspec = spectrum * vmask
  bdrspec *= vmask
  
  ;; For considering the HISA spectrum (as opposed to the overall HI
  ;; spectrum), we mask in regions w/ bdrspec > (TA = !MW.TA_HI)
  ;; (Somewhat following Gibson et al., 2005, ApJ, 626,214 -- they use
  ;; TA = 70k)).
  ;; We use the total region enclosed by [first,last] time
  ;; spectrum goes above this level.  This excludes fluctuations in
  ;; the thin emission at the extremes of the velocity spectrum  from
  ;; being interpreted as significant HISA.
  maskind = WHERE( bdrspec GE !MW.TA_HI, nmask )
  emmask = bdrspec * 0
  IF nmask NE 0 THEN $
     emmask[maskind[0]:maskind[nmask-1]] = 1
  
  ;; Look only for where absorption is present
  spectrum = ((spectrum - bdrspec)<0.) * emmask
  
  ;; Deal with NaN's...
  ind = WHERE( ~ FINITE(spectrum), nind)
  IF nind NE 0 THEN spectrum[ind] = 0
  
  ;;====================================================================
  ;; Plot, if desired
  IF KEYWORD_SET( plot ) THEN BEGIN
     plot,v_std,onspec,/xst,xtitle='V!DLSR!N (km/s)',$
          ytitle='!6T!DA!N (K)',linestyle=2, yrange=[0,180],/yst,$
          title='BGPS #'+string(bgps.cnum,format="(I04)")
     oplot,v_std,bdrspec,ps=10,color='0000ff'x,linestyle=2
     
     ;; Oplot the masks
     oplot,v_std,emmask * !MW.TA_HI,ps=10,color='00ffff'x,linestyle=3
     oplot,v_std,vmask*100,ps=10,color='ff8915'x,linestyle=4
     
     ;; Plot the HISA spectrum in a visible manner...
     oplot,v_std,spectrum + 170 ,ps=10,color='ffff00'x
     
     ;; Velocity limits (to exclude significant noise at large V_LSR).
     vline,v_max
     vline,v_max+20.,color='80ff80'x
     al_legend,position=[-90,150],['onspec','bdrspec'],linestyle=2,$
               color=['ffffff'x,'0000ff'x]
     
  ENDIF
  ;;====================================================================

  ;; In Erik's GRSMATCH.pro, he squared the output spectrum
  ;; just before returning it to the calling routine.  Not sure WHY he
  ;; did this, other than to make peakier the spectrum weights for the
  ;; kinematic distance pdfs.
  ;; Changing 4/29/11, HISA.pro is no longer squaring the output
  ;; spectrum.  Since we enforce looking only for ABSOPRTION in the
  ;; spectrum, we can take the absolute value of the output spectrum
  ;; to yield our requisite spectrum.
  ;; 
  ;; We wish to NOT square the spectrum so that we can directly use
  ;; the antenna temperature decrement for other purposes, and so that
  ;; the WIDTHS of the spectral features are physical and directly
  ;; comparable with the dense-gas spectra.

  spectrum = ABS(spectrum)
  
  ;; If blank, simply return null_spec
  IF total(spectrum) EQ 0 THEN RETURN,null_spec
  
  RETURN, spectrum
END
